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We extend the treatment of the problem of the gravitational waves produced by a particle of 
negligible mass orbiting a Kerr black hole using black hole perturbation theory [17] in the time 
domain, to elliptic and inclined orbits. We model the particle by smearing the singularities in 
the source term using narrow Gaussian distributions. We compare results (energy and angular 
momentum fluxes) for such orbits with those computed using the frequency domain formalism. 
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I. INTRODUCTION 

Current gravitational wave detectors such as LIGO, VIRGO, GEO and TAMA are now starting "science runs" and 
the plans for upgrades to these detectors are already in the works. This paper attempts to model a likely occurrence 
in galactic nuclei i.e. the capture of compact stellar-sized objects by the supermassive black holes that exist in the 
center of most galaxies. These type of events will generate gravitational waves within the sensitivity band of NASA's 
space-based LISA project. 

This problem has been treated in the past, quite extensively, using a frequency domain approach to linear black 
hole perturbation theory [1-16]. A similar calculation in the time domain was only recently presented for the special 
case of equatorial, circular particle orbits around a rapidly rotating black hole [17]. In this work, we extend the time 
domain calculation to treat more generic types of orbits, i.e. elliptic and inclined, and compare energy and angular 
momentum fluxes lost from gravitational wave emission to current frequency domain calculations. We ignore the 
effects of radiation reaction in this work, with the hope of addressing those in detail sometime in the future. We note 
that for highly elliptic orbits, both approaches develop problems and only a very detailed study, possibly with the 
inclusion of radiation reaction effects will lead us to favor one over the other. 

This paper is organized as follows. In the next section we review, briefly, the basics of our time domain approach. 
Details can be found in the second section of our past work [17]. Then in the following section, we present our results 
for waveforms and radiative fluxes and compare them with frequency domain results. In the last section, we present 
our conclusions and emphasize that future work will be based on the inclusion of radiation reaction in our calculations. 

II. TEUKOLSKY EQUATION IN THE TIME DOMAIN 

The general approach to this problem in the time domain, with details regarding numerical implementation, stability, 
second-order convergence, etc. have been presented in the earlier work on the subject [17], so we simply summarize 
some of the essentials here. The Teukolsky equation in Boyer-Lindquist coordinates with the matter-source term, 
appears below [18], 
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where A = 

The methodology for numerically evolving the homogeneous part of this equation has been presented in detail at 
various places in the literature [19]. In this work we focus on the matter-source term. In the time domain Teukolsky 
equation this source term is given by 
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All the greek symbols used above are the ones as defined in the original Teukolsky equation paper [18]. By choosing 
their values in Boyer-Lindquist coordinates and expanding everything we obtain an explicit but quite complicated 
form of the source term. We treat the delta functions that appear in the energy-momentum tensor of the particle in 
orbit by using an approximation in which we substitute the delta function by a very narrow Gaussian function (whose 
width is only of a few grid points), 
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We handle the deltas in the radial and angular direction through the Gaussian approximation, whereas the d{<p — (p(t)) 
function we handle analytically. We use code generation software to convert the source term expressions to FORTRAN 
code that can be used with our Teukolsky equation evolver. We verify that this substitution is adequate by showing 
convergence of the obtained waveform and radiative flux values as cr ^ 0. We also take special care that the Gaussian 
integrates to unity on the spatial hypersurface. To achieve this with the three-metric of the background black hole 
spacetime, we normalize the mass-density of the particle by the factor. 
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where 7^^^ is the 3-metric of flat space and g^^^ is the 3-metric of a slice of constant Boyer-Lindquist time of the Kerr 
spacetime. 

The initial data we take for the fields is zero. Since we are considering the Teukolsky equation with a source this 
would correspond to the particle "appearing suddenly" and that generates an artificial burst of radiation. We handle 
this by evolving to late times and computing fluxes only when the system has settled down. Sometime in the future, 
we will examine the implementation of initial data in this context, but at this time we will simply rid the system of 
this initial burst, by evolving "long enough". 

Lastly, we also need a high accuracy general geodesic equation integrator, to obtain the position of the particle at 
every time step. The main problem with integrating the geodesic equations is the presence of turning points. At the 
turning points, since the right-hand sides of the geodesic equations go to zero (see below), the integrators have trouble 
resolving the motion correctly and large amounts of error accumulate. 
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Fortunately, there is such a public domain integrator, called "Good" written by Scott Hughes [20]. In Geod, this 
problem is solved by transforming the equations to orbital parameters that do not have turning points. More in- 
formation can be found about this approach in the documentation of this software package. We call this geodesic 
equation integrator from within our Teukolsky code to update the particle's position at every time step. 



III. RESULTS 



In this section, we present sample waveforms and compute radiative fluxes from numerical evolutions of the Teukol- 
sky equation with matter-source term as discussed in the previous section. As detailed before [17], the numerical 

implementation of this work is in the form of a set of 2+1 dimensional linear PDE's that uses the Lax-Wendroff 
evolution scheme. Also discussed therein is the convergence and stability of the entire approach, therefore that study 
shall not be presented here. 
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m mode 


Time domain energy flux 


Frequency domain flux 


1 


2.9 X 10"'" 


2.8 X 10"'" 


2 


1.4 X 10""' 


1.8 X 10""' 


3 


4.2 X IQ-"" 


5.3 X IQ-"" 



TABLE I: Comparisons of gravitational wave energy fluxes detected at infinity using the frequency domain solution with the 
results measured mimcrically at r = lOOM while evolving the Teukolsky equation in the time domain as discussed here. The 
particle here is in an elliptic, equatorial orbit of semi-latus rectum 2.0 rjsco and eccentricity of 0.5 about a black hole with Kerr 
parameter a/M = 0.9 and n/M = 0.01. 



We calculate the energy and angular momentum fluxes emitted due to gravity wave emission from the binary for 
a sample equatorial, elliptic orbit and show close agreement with current results of high accuracy frequency domain 

calculations [22]. The code monitors the value of the Teukolsky fTinction at a set radial location in the grid far from the 
horizon and computes the energy and angular momentum fluxes in gravitational radiation according to the formulae 
[21]: 
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In Tables I and II we present representative comparisons with frequency domain calculations for equatorial, elliptic 
orbits. The level of agreement between the two approaches is about the same (within 25%) over a wide range 
of eccentricities (0.0 through 0.8). It is interesting to note both frequency domain approaches and time domain 
approaches begin to have computational issues for orbits of high eccentricity (above 0.8) and also for hyperbolic 
orbits. In the frequency domain, where the method is based on variable separating the Teukolsky equation the need 
to compute for a large number of k and i modes for a given accuracy, makes it computationally very demanding 
(see Figure 1). In the time domain case, since highly elliptic orbits have large orbital periods, we need to evolve for 
very long times (several tens of thousands times the mass of the hole) in order to compute average radiation fluxes 
(see Figure 3). Therefore it is not clear at this point, which of the two approaches is better suited to address the 
general problem of extreme mass ratio inspiral. However, this question will be revisited, once more work is done on 
the inclusion of radiation reaction and a more detailed comparison study is undertaken. 

We also include a representative waveform in Figure 2 from the evolution of a particle in an equatorial, elliptic orbit. 
It shows the real part of the m = 3 mode of the Teukolsky function, and as mentioned before, the initial burst of 
radiation is unphysical, since it corresponds to the particle "appearing out of nowhere" . It is only after the evolution 
has settled down to a time-periodic pattern, do we compute average fliixcs. In Figure 3 we show the variation of 
energy fluxes with time, from the evolution of a particle in an equatorial, elliptic orbit and also in an inclined, circular 
orbit. Its is interesting to note that if we keep the semi-latus rectum fixed, and study the energy-fluxes as they change 
with eccentricity for elliptic-equatorial orbits, they increase at first, but then begin to drop after reaching a maximum 
at an eccentricity of about 0.5. It turns out that this variation can be understood easily by studying the same in the 
context of a binary Newtownian system where this turn-over point can be computed to be about 0.47 [23]. 

Lastly, we suggest here an alternate approximate approach to treating delta functions on a spatial grid. This 
approach is significantly more efficient (more than two orders of magnitude faster in a numerical implementation) 
and yields significantly more accurate results (in terms of their agreement with frequency domain) compared with 
the above mentioned narrow Gaussian approach. One can think of this as a limiting case of the Gaussian particle 
approach, where the width of the Gaussian is made so narrow that it is comparable to the numerical grid separation. 
In this approach, we implement a "discrete" version of the delta function on the spatial grid. More specifically, we 
define a discrete-delta function as being unity at the appropriate grid point and zero everywhere else. We then define 
the derivatives of the delta function as the result of applying the usual finite-difference derivative stencils to the above 
mentioned discrete-delta fimction. Of course, the discrete-delta needs to be normalized suitably, i.e. scaled by a factor 
such that its integral on the spatial grid is unity. Note that such a delta function satisfies the usual delta function 
properties, such as f{x)6'{x) = —f'{x)S{x), on the discrete grid. Care has to be taken to study convergence of the 
code using such an approach, however, just the efficiency and accuracy obtained by using this simple idea, makes it 
worth mentioning and investigating further. 
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m mode 


Time domain ang. mom. flux 


Frequency domain flux 


1 


2.4 X 10""" 


2.3 X 10"™ 


2 


1.0 X 10""" 


1.4 X 10""" 


3 


3.0 X 10""' 


4.0 X 10""' 



TABLE II: Comparisons of gravitational wave angular momentum fluxes detected at infinity using the frequency domain 
solution with the results measured numerically at r = 100i\f while evolving the Teukolsky equation in the time domain as 
discussed here. The particle here is in a elliptic, equatorial orbit of semi-latus rectum 2.0 Visco and eccentricity of 0.5 about a 
black hole with Kerr parameter a/M = 0.9 and i-i/M = 0.01. 
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FIG. 1: Energy fluxes computed from the evolution of the m — 3 mode of the Teukolsky function from a particle (i-i/M — 0.01) 
orbiting a rapidly rotating Kerr hole of a = 0.9M. These results were obtained from a frequency domain computation, and 
are plotted against the k modes. All the orbits are equatorial and have a semi-latus rectum of 2ri3co- Note that as the orbit's 
eccentricity increases, results from more k modes are necessary, increasing the computational needs of the problem. 



IV. CONCLUSIONS 



We can now treat different types of particle orbits in the context of extreme mass ratio inspiral, using black hole 
perturbation theory in the time domain. We get good agreement in comparisons with current frequency domain 
results for energy and angular momentum radiative fluxes. At this time, it is not clear which of these two approaches 
is better suited to the generic case and will provide more accurate results, more efficiently. 

In future work, we hope to be able to further widen the scope of these calculations by incorporating the case of a 
spinning particle in orbit about a Kerr hole and we also expect to investigate methods to include radiation reaction. 
The first approach will be to correct the orbit trajectory using energy-balance arguments, as for instance in [14]. But, 
later other proposals for the inclusion of radiation reaction forces will be incorporated. 
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FIG. 2: Representative data showing the evolution of the real part of the Tcukolsky function for the m — 3 mode of a particle 
{n/M = 0.01) orbiting a rapidly rotating Kerr hole of a = 0.9M in an elliptic orbit of semi-latus rectum 2 risco {risco = 2.32M) 
and eccentricity of 0.5 at an observation point r/M = 100, 9 = ■k/2. 
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FIG. 3: Energy fluxes computed from the evolution of the m = 3 mode of the Teukolsky function from a particle {n/M = 0.01) 
orbiting a rapidly rotating Kerr hole of a = O.QAf. The top plot is from the particle in an equatorial, elliptic orbit of somi-latus 
rectum 2risco {risco ~ 2.32M) and eccentricity of 0.5 at an observation point r/M = 100. It should be clear from this plot, 
that as eccentricity increases, longer evolutions are needed to compute average fluxes. The bottom plot is from the particle in 
a circular, inclined orbit with radius 2 risco and initial inclination angle of 7r/4 radians. 



